home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / betabase.c < prev    next >
C/C++ Source or Header  |  1990-10-04  |  8KB  |  300 lines
  1. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  2. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  3. /* You may give out copies of this software; for conditions see the    */
  4. /* file COPYING included with this distribution.                       */
  5.  
  6. #include "xlisp.h"
  7. #include "osdef.h"
  8. #ifdef ANSI
  9. #include "xlsproto.h"
  10. #else
  11. #include "xlsfun.h"
  12. #endif ANSI
  13.  
  14. /* forward declarations */
  15. #ifdef ANSI
  16. double logbeta(double,double),betai(double,double,double),
  17.        xinbta(double *,double *,double *,double *,int *);
  18. #else
  19. double logbeta(),betai(),xinbta();
  20. #endif ANSI
  21.  
  22. void betabase(x, a, b, gia, gib, cdf)
  23.      int *gia, *gib;
  24.      double *x, *a, *b, *cdf;
  25. {
  26.   *cdf = betai(*x, *a, *b);
  27. }
  28.  
  29. double ppbeta(p, a, b, ifault)
  30.      double p, a, b;
  31.      int *ifault;
  32. {
  33.   double lbeta;
  34.  
  35.   lbeta = gamma(a) + gamma(b) - gamma(a + b);
  36.   return(xinbta(&a, &b, &lbeta, &p, ifault));
  37. }
  38.  
  39. /*
  40.   Static routines
  41. */
  42.  
  43. static double logbeta(p, q)
  44.      double p, q;
  45. {
  46.   return(gamma(p) + gamma(q) - gamma(p + q));
  47. }
  48. /* defined in headers JKL
  49. #define Min(x,y) (((x) < (y)) ? (x) : (y))
  50. #define Max(x,y) (((x) > (y)) ? (x) : (y))
  51. */
  52. static double betai(x, pin, qin)
  53.      double x, pin, qin;
  54. {
  55.   /* Translated from FORTRAN
  56.      july 1977 edition.  w. fullerton, c3, los alamos scientific lab.
  57.      based on bosten and battiste, remark on algorithm 179, comm. acm,
  58.      v 17, p 153, (1974).
  59.      
  60.      input arguments --
  61.      x      upper limit of integration.  x must be in (0,1) inclusive.
  62.      p      first beta distribution parameter.  p must be gt 0.0.
  63.      q      second beta distribution parameter.  q must be gt 0.0.
  64.      betai  the incomplete beta function ratio is the probability that a
  65.             random variable from a beta distribution having parameters
  66.             p and q will be less than or equal to x.
  67.   */
  68.   double c, finsum, p, ps, q, term, xb, xi, y, dbetai, p1;
  69.   static double eps = 0.0, alneps = 0.0, sml = 0.0, alnsml = 0.0;
  70.   int i, n, ib;
  71.  
  72.   /* I'm not sure these tolerances are optimal */
  73.   if (eps == 0.0) {
  74.     eps = macheps();
  75.     alneps = log(eps);
  76.     sml = macheps();
  77.     alnsml = log(sml);
  78.   }
  79.  
  80.   y = x;
  81.   p = pin;
  82.   q = qin;
  83.   if (q > p || x >= 0.8)
  84.     if (x >= 0.2) {
  85.       y = 1.0 - y;
  86.       p = qin;
  87.       q = pin;
  88.     }
  89.  
  90.   if ((p + q) * y / (p + 1.0) < eps) {
  91.     dbetai = 0.0;
  92.     xb = p * log(Max(y, sml)) - log(p) - logbeta(p, q);
  93.     if (xb > alnsml && y != 0.0) dbetai = exp(xb);
  94.     if (y != x || p != pin) dbetai = 1.0 - dbetai;
  95.   }
  96.   else {
  97.     /*
  98.      *  evaluate the infinite sum first.  term will equal
  99.      *  y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
  100.      */
  101.     ps = q - floor(q);
  102.     if (ps == 0.0) ps = 1.0;
  103.     xb = p * log(y) - logbeta(ps, p) - log(p);
  104.     dbetai = 0.0;
  105.     if (xb >= alnsml) {
  106.  
  107.       dbetai = exp(xb);
  108.       term = dbetai * p;
  109.       if (ps != 1.0) {
  110.         n = Max(alneps / log(y), 4.0);
  111.         for (i = 1; i <= n; i++) {
  112.           xi = i;
  113.           term = term * (xi - ps) * y / xi;
  114.           dbetai = dbetai + term / (p + xi);
  115.         }
  116.       }
  117.     }
  118.     /*
  119.      * now evaluate the finite sum, maybe.
  120.      */
  121.     if (q > 1.0) {
  122.  
  123.       xb = p * log(y) + q * log(1.0 - y) - logbeta(p,q) - log(q);
  124.       ib = Max(xb / alnsml, 0.0);
  125.       term = exp(xb - ((float) ib) * alnsml);
  126.       c = 1.0 / (1.0 - y);
  127.       p1 = q * c / (p + q - 1.0);
  128.  
  129.       finsum = 0.0;
  130.       n = q;
  131.       if (q == (float) n) n = n - 1;
  132.       for (i = 1; i <= n; i++) {
  133.         if (p1 <= 1.0 && term / eps <= finsum) break;
  134.         xi = i;
  135.         term = (q - xi + 1.0) * c * term / (p + q - xi);
  136.  
  137.         if (term > 1.0) ib = ib - 1;
  138.         if (term > 1.0)    term = term * sml;
  139.  
  140.         if (ib==0) finsum = finsum + term;
  141.       }
  142.  
  143.       dbetai = dbetai + finsum;
  144.     }
  145.     if (y != x || p != pin) dbetai = 1.0 - dbetai;
  146.     dbetai = Max(Min(dbetai, 1.0), 0.0);
  147.   }
  148.   return(dbetai);
  149. }
  150.  
  151. /*
  152.   xinbta.f -- translated by f2c and modified
  153.   
  154.   algorithm as 109 appl. statist. (1977), vol.26, no.1
  155.   (replacing algorithm as 64  appl. statist. (1973), vol.22, no.3)
  156.  
  157.   Remark AS R83 has been incorporated in this version.
  158.  
  159.   Computes inverse of the incomplete beta function
  160.   ratio for given positive values of the arguments
  161.   p and q, alpha between zero and one.
  162.   log of complete beta function, beta, is assumed to be known.
  163.  
  164.   Auxiliary function required: betai
  165.  
  166.   SAE below is the most negative decimal exponent which does not
  167.   cause an underflow; a value of -308 or thereabouts will often be
  168. */
  169.  
  170. static double xinbta(p, q, beta, alpha, ifault)
  171.      double *p, *q, *beta, *alpha;
  172.      int *ifault;
  173. {
  174.   /* Initialized data */
  175.   static double sae = -30.0; /* this should be sufficient */
  176.   static double zero = 0.0;
  177.   static double one = 1.0;
  178.   static double two = 2.0;
  179.   static double three = 3.0;
  180.   static double four = 4.0;
  181.   static double five = 5.0;
  182.   static double six = 6.0;
  183.  
  184.   /* System generated locals */
  185.   double ret_val, d_1, d_2;
  186.  
  187.   /* Local variables */
  188.   static int indx;
  189.   static double prev, a, g, h, r, s, t, w, y, yprev, pp, qq;
  190.   static double sq, tx, adj, acu;
  191.   static int iex;
  192.   static double fpu, xin;
  193.  
  194.   /* Define accuracy and initialise. */
  195.   fpu = sae * 10.;
  196.   ret_val = *alpha;
  197.  
  198.   /* test for admissibility of parameters */
  199.   *ifault = 1;
  200.   if (*p <= zero || *q <= zero) return ret_val;
  201.   *ifault = 2;
  202.   if (*alpha < zero || *alpha > one) return ret_val;
  203.   *ifault = 0;
  204.   if (*alpha == zero || *alpha == one) return ret_val;
  205.  
  206.   /* change tail if necessary */
  207.   if (*alpha <= .5) {
  208.     a = *alpha;
  209.     pp = *p;
  210.     qq = *q;
  211.     indx = FALSE;
  212.   }
  213.   else {
  214.     a = one - *alpha;
  215.     pp = *q;
  216.     qq = *p;
  217.     indx = TRUE;
  218.   }
  219.  
  220.   /* calculate the initial approximation */
  221.   r = sqrt(-log(a * a));
  222.   y = r - (r * .27061 + 2.30753) / (one + (r * .04481 + .99229) * r);
  223.   if (pp > one && qq > one) {
  224.     r = (y * y - three) / six;
  225.     s = one / (pp + pp - one);
  226.     t = one / (qq + qq - one);
  227.     h = two / (s + t);
  228.     d_1 = y * sqrt(h + r) / h;
  229.     d_2 = (t - s) * (r + five / six - two / (three * h));
  230.     w = d_1 - d_2;
  231.     ret_val = pp / (pp + qq * exp(w + w));
  232.   }
  233.   else {
  234.     r = qq + qq;
  235.     t = one / (qq * 9.);
  236.     /* Computing 3rd power */
  237.     d_1 = one - t + y * sqrt(t);
  238.     t = r * (d_1 * d_1 * d_1);
  239.     if (t <= zero) {
  240.       ret_val = one - exp((log((one - a) * qq) + *beta) / qq);
  241.     }
  242.     else {
  243.       t = (four * pp + r - two) / t;
  244.       if (t <= one) ret_val = exp((log(a * pp) + *beta) / pp);
  245.       else ret_val = one - two / (t + one);
  246.     }
  247.   }
  248.  
  249.  
  250.   /* 
  251.     solve for x by a modified newton-raphson method, using the function betai
  252.   */
  253.   r = one - pp;
  254.   t = one - qq;
  255.   yprev = zero;
  256.   sq = one;
  257.   prev = one;
  258.   if (ret_val < 1e-4) ret_val = 1e-4;
  259.   if (ret_val > .9999) ret_val = .9999;
  260.   /* Computing MAX, two 2nd powers */
  261.   d_1 = -5.0 / (pp * pp) - 1.0 / (a * a) - 13.0;
  262.   iex = (sae > d_1) ? sae : d_1;
  263.   acu = pow(10.0, (double) iex);
  264.   do {
  265.     y = betai(ret_val, pp, qq);
  266.     if (*ifault != 0) {
  267.       *ifault = 3;
  268.       return ret_val;
  269.     }
  270.     xin = ret_val;
  271.     y = (y - a) * exp(*beta + r * log(xin) + t * log(one - xin));
  272.     if (y * yprev <= zero) {
  273.       prev = (sq > fpu) ? sq : fpu;
  274.     }
  275.     g = one;
  276.     do {
  277.       adj = g * y;
  278.       sq = adj * adj;
  279.       if (sq < prev) {
  280.         tx = ret_val - adj;
  281.         if (tx >= zero && tx <= one) {
  282.           if (prev <= acu || y * y <= acu) {
  283.             if (indx) ret_val = one - ret_val;
  284.             return ret_val;
  285.           }
  286.           if (tx != zero && tx != one) break;
  287.         }
  288.       }
  289.       g /= three;
  290.     } while (TRUE);
  291.     if (tx == ret_val) {
  292.       if (indx) ret_val = one - ret_val;
  293.       return ret_val;
  294.     }
  295.     ret_val = tx;
  296.     yprev = y;
  297.   } while (TRUE);
  298.   return ret_val;
  299. }
  300.